In [ ]:
import stim
import pymatching
import numpy as np
import sinter
from typing import List
import matplotlib.pyplot as plt
In [ ]:
def generate_circuit(rounds=1, before_round_data_depolarization= 0.1, before_measure_flip_probability=0.2, after_reset_flip_probability=0.3, after_clifford_depolarization=0.4):
    return f"""# surface_code circuit exercise.
# rounds: {rounds}
# before_round_data_depolarization: {before_round_data_depolarization}
# before_measure_flip_probability: {before_measure_flip_probability}
# after_reset_flip_probability: {after_reset_flip_probability}
# after_clifford_depolarization: {after_clifford_depolarization}
# layout:
# L00 X01 L02 X03 L04 X05 L06
# Z10 d11 Z12 d13 Z14 d15 Z16
# d20 X21 d22 X22 d24 X25 d26
# Z30 d31 Z32 d33 Z34 d35 Z36
# d40 X41 d42 X42 d44 X45 d46
# Z50 d51 Z52 d53 Z54 d55 Z56
# d60 X61 d62 X62 d64 X65 d66
# Legend:
#     d# = data qubit
#     L# = data qubit with logical observable crossing
#     X# = measurement qubit (X stabilizer)
#     Z# = measurement qubit (Z stabilizer)
QUBIT_COORDS(0, 0) 1
QUBIT_COORDS(0, 1) 2
QUBIT_COORDS(0, 2) 3
QUBIT_COORDS(0, 3) 4
QUBIT_COORDS(0, 4) 5
QUBIT_COORDS(0, 5) 6
QUBIT_COORDS(0, 6) 7 
QUBIT_COORDS(1, 0) 8
QUBIT_COORDS(1, 1) 9
QUBIT_COORDS(1, 2) 10
QUBIT_COORDS(1, 3) 11
QUBIT_COORDS(1, 4) 12
QUBIT_COORDS(1, 5) 13
QUBIT_COORDS(1, 6) 14
QUBIT_COORDS(2, 0) 15
QUBIT_COORDS(2, 1) 16
QUBIT_COORDS(2, 2) 17
QUBIT_COORDS(2, 3) 18
QUBIT_COORDS(2, 4) 19
QUBIT_COORDS(2, 5) 20
QUBIT_COORDS(2, 6) 21
QUBIT_COORDS(3, 0) 22
QUBIT_COORDS(3, 1) 23
QUBIT_COORDS(3, 2) 24
QUBIT_COORDS(3, 3) 25
QUBIT_COORDS(3, 4) 26
QUBIT_COORDS(3, 5) 27
QUBIT_COORDS(3, 6) 28
QUBIT_COORDS(4, 0) 29
QUBIT_COORDS(4, 1) 30
QUBIT_COORDS(4, 2) 31
QUBIT_COORDS(4, 3) 32
QUBIT_COORDS(4, 4) 33
QUBIT_COORDS(4, 5) 34
QUBIT_COORDS(4, 6) 35
QUBIT_COORDS(5, 0) 36
QUBIT_COORDS(5, 1) 37
QUBIT_COORDS(5, 2) 38
QUBIT_COORDS(5, 3) 39
QUBIT_COORDS(5, 4) 40
QUBIT_COORDS(5, 5) 41
QUBIT_COORDS(5, 6) 42
QUBIT_COORDS(6, 0) 43
QUBIT_COORDS(6, 1) 44
QUBIT_COORDS(6, 2) 45
QUBIT_COORDS(6, 3) 46
QUBIT_COORDS(6, 4) 47
QUBIT_COORDS(6, 5) 48
QUBIT_COORDS(6, 6) 49
R 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
X_ERROR({after_reset_flip_probability}) 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
R 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
X_ERROR({after_reset_flip_probability}) 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
TICK
DEPOLARIZE1({before_round_data_depolarization}) 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
TICK
# Plaquette Syndromes 
H 2 4 6 16 18 20 30 32 34 44 46 48
DEPOLARIZE1({after_clifford_depolarization}) 2 4 6 16 18 20 30 32 34 44 46 48
TICK
# Plaquette Syndromes triplets
CX 1 2 3 4 5 6 43 44 45 46 47 48
DEPOLARIZE2({after_clifford_depolarization}) 1 2 3 4 5 6 43 44 45 46 47 48
TICK
CX 9 2 11 4 13 6 37 44 39 46 41 48
DEPOLARIZE2({after_clifford_depolarization}) 9 2 11 4 13 6 37 44 39 46 41 48
TICK
CX 3 2 5 4 7 6 45 44 47 46 49 48
DEPOLARIZE2({after_clifford_depolarization}) 3 2 5 4 7 6 45 44 47 46 49 48
TICK
# Plaquette Syndromes quads
CX 9 16 11 18 13 20 23 30 25 32 27 34
DEPOLARIZE2({after_clifford_depolarization}) 9 16 11 18 13 20 23 30 25 32 27 34
TICK
CX 15 16 17 18 19 20 29 30 31 32 33 34
DEPOLARIZE2({after_clifford_depolarization}) 15 16 17 18 19 20 29 30 31 32 33 34
TICK
CX 17 16 19 18 21 20 31 30 33 32 35 34
DEPOLARIZE2({after_clifford_depolarization}) 17 16 19 18 21 20 31 30 33 32 35 34
TICK
CX 23 16 25 18 27 20 37 30 39 32 41 34
DEPOLARIZE2({after_clifford_depolarization}) 23 16 25 18 34 20 37 30 39 32 41 34
TICK
# Vertex Syndromes triplets
CX 1 8 15 22 29 36 7 14 21 28 35 42
DEPOLARIZE2({after_clifford_depolarization}) 1 8 15 22 29 36 7 14 21 28 35 42
TICK 
CX 15 8 29 22 43 36 21 14 35 28 49 42
DEPOLARIZE2({after_clifford_depolarization}) 1 8 15 22 29 36 7 14 21 28 35 42
TICK 
CX 9 8 23 22 37 36 13 14 27 28 41 42
DEPOLARIZE2({after_clifford_depolarization}) 1 8 15 22 29 36 7 14 21 28 35 42
TICK 
# Vertex Syndromes quads
CX 9 10 11 12 23 24 25 26 37 38 39 40
DEPOLARIZE2({after_clifford_depolarization}) 9 10 11 12 23 24 25 26 37 38 39 40
TICK
CX 3 10 5 12 17 24 19 26 31 38 33 40
DEPOLARIZE2({after_clifford_depolarization}) 3 10 5 12 17 24 19 26 31 38 33 40
TICK
CX 11 10 13 12 25 24 27 26 39 38 41 40
DEPOLARIZE2({after_clifford_depolarization}) 11 10 13 12 25 24 27 26 39 38 41 40
TICK
CX 17 10 19 12 31 24 33 26 45 38 47 40
DEPOLARIZE2({after_clifford_depolarization}) 17 10 19 12 31 24 33 26 45 38 47 40
TICK
# Plaquette Syndromes 
H 2 4 6 16 18 20 30 32 34 44 46 48
DEPOLARIZE1({after_clifford_depolarization}) 2 4 6 16 18 20 30 32 34 44 46 48
TICK
# Measure measurement qubits
X_ERROR({before_measure_flip_probability}) 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
MR 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
X_ERROR({after_reset_flip_probability}) 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
# Detect errors using the vertex Syndromes
DETECTOR(1, 0, 0) rec[-21]
DETECTOR(3, 0, 0) rec[-14]
DETECTOR(5, 0, 0) rec[-7]
DETECTOR(1, 2, 0) rec[-20]
DETECTOR(3, 2, 0) rec[-13]
DETECTOR(5, 2, 0) rec[-6]
DETECTOR(1, 4, 0) rec[-19]
DETECTOR(3, 4, 0) rec[-12]
DETECTOR(3, 4, 0) rec[-5]
DETECTOR(1, 6, 0) rec[-18]
DETECTOR(3, 6, 0) rec[-11]
DETECTOR(5, 6, 0) rec[-4]
REPEAT {rounds}""" + "{" + f"""
    TICK
    DEPOLARIZE1({before_round_data_depolarization}) 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
    TICK
    # Plaquette Syndromes 
    H 2 4 6 16 18 20 30 32 34 44 46 48
    DEPOLARIZE1({after_clifford_depolarization}) 2 4 6 16 18 20 30 32 34 44 46 48
    TICK
    # Plaquette Syndromes triplets
    CX 1 2 3 4 5 6 43 44 45 46 47 48
    DEPOLARIZE2({after_clifford_depolarization}) 1 2 3 4 5 6 43 44 45 46 47 48
    TICK
    CX 9 2 11 4 13 6 37 44 39 46 41 48
    DEPOLARIZE2({after_clifford_depolarization}) 9 2 11 4 13 6 37 44 39 46 41 48
    TICK
    CX 3 2 5 4 7 6 45 44 47 46 49 48
    DEPOLARIZE2({after_clifford_depolarization}) 3 2 5 4 7 6 45 44 47 46 49 48
    TICK
    # Plaquette Syndromes quads
    CX 9 16 11 18 13 20 23 30 25 32 27 34
    DEPOLARIZE2({after_clifford_depolarization}) 9 16 11 18 13 20 23 30 25 32 27 34
    TICK
    CX 15 16 17 18 19 20 29 30 31 32 33 34
    DEPOLARIZE2({after_clifford_depolarization}) 15 16 17 18 19 20 29 30 31 32 33 34
    TICK
    CX 17 16 19 18 21 20 31 30 33 32 35 34
    DEPOLARIZE2({after_clifford_depolarization}) 17 16 19 18 21 20 31 30 33 32 35 34
    TICK
    CX 23 16 25 18 27 20 37 30 39 32 41 34
    DEPOLARIZE2({after_clifford_depolarization}) 23 16 25 18 34 20 37 30 39 32 41 34
    TICK
    # Vertex Syndromes triplets
    CX 1 8 15 22 29 36 7 14 21 28 35 42
    DEPOLARIZE2({after_clifford_depolarization}) 1 8 15 22 29 36 7 14 21 28 35 42
    TICK 
    CX 15 8 29 22 43 36 21 14 35 28 49 42
    DEPOLARIZE2({after_clifford_depolarization}) 1 8 15 22 29 36 7 14 21 28 35 42
    TICK 
    CX 9 8 23 22 37 36 13 14 27 28 41 42
    DEPOLARIZE2({after_clifford_depolarization}) 1 8 15 22 29 36 7 14 21 28 35 42
    TICK 
    # Vertex Syndromes quads
    CX 9 10 11 12 23 24 25 26 37 38 39 40
    DEPOLARIZE2({after_clifford_depolarization}) 9 10 11 12 23 24 25 26 37 38 39 40
    TICK
    CX 3 10 5 12 17 24 19 26 31 38 33 40
    DEPOLARIZE2({after_clifford_depolarization}) 3 10 5 12 17 24 19 26 31 38 33 40
    TICK
    CX 11 10 13 12 25 24 27 26 39 38 41 40
    DEPOLARIZE2({after_clifford_depolarization}) 11 10 13 12 25 24 27 26 39 38 41 40
    TICK
    CX 17 10 19 12 31 24 33 26 45 38 47 40
    DEPOLARIZE2({after_clifford_depolarization}) 17 10 19 12 31 24 33 26 45 38 47 40
    TICK
    # Plaquette Syndromes 
    H 2 4 6 16 18 20 30 32 34 44 46 48
    DEPOLARIZE1({after_clifford_depolarization}) 2 4 6 16 18 20 30 32 34 44 46 48
    TICK
    # Measure measurement qubits
    X_ERROR({before_measure_flip_probability}) 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
    MR 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
    X_ERROR({after_reset_flip_probability}) 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
    # Used for adding time?
    SHIFT_COORDS(0, 0, 1)
    # Detecting errors by comparing to previous shot
    DETECTOR(6, 5, 0) rec[-1] rec[-25]
    DETECTOR(6, 3, 0) rec[-2] rec[-26]
    DETECTOR(6, 1, 0) rec[-3] rec[-27]
    DETECTOR(5, 6, 0) rec[-4] rec[-28]
    DETECTOR(5, 4, 0) rec[-5] rec[-29]
    DETECTOR(5, 2, 0) rec[-6] rec[-30]
    DETECTOR(5, 0, 0) rec[-7] rec[-31]
    DETECTOR(4, 5, 0) rec[-8] rec[-32]
    DETECTOR(4, 3, 0) rec[-9] rec[-33]
    DETECTOR(4, 1, 0) rec[-10] rec[-34]
    DETECTOR(3, 6, 0) rec[-11] rec[-35]
    DETECTOR(3, 4, 0) rec[-12] rec[-36]
    DETECTOR(3, 2, 0) rec[-13] rec[-37]
    DETECTOR(3, 0, 0) rec[-14] rec[-38]
    DETECTOR(2, 5, 0) rec[-15] rec[-39]
    DETECTOR(2, 3, 0) rec[-16] rec[-40]
    DETECTOR(2, 1, 0) rec[-17] rec[-41]
    DETECTOR(1, 6, 0) rec[-18] rec[-42]
    DETECTOR(1, 4, 0) rec[-19] rec[-43]
    DETECTOR(1, 2, 0) rec[-20] rec[-44]
    DETECTOR(1, 0, 0) rec[-21] rec[-45]
    DETECTOR(0, 5, 0) rec[-22] rec[-46]
    DETECTOR(0, 3, 0) rec[-23] rec[-47]
    DETECTOR(0, 1, 0) rec[-24] rec[-48]
"""+"}"+f"""
# Final readout of the 25 data qubits in Z basis as we are using vertex syndromes
X_ERROR({before_measure_flip_probability}) 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
M 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
DETECTOR(1, 0, 1) rec[-18] rec[-21] rec[-25] rec[-46]
DETECTOR(3, 0, 1) rec[-11] rec[-14] rec[-18] rec[-39]
DETECTOR(5, 0, 1) rec[-4] rec[-7] rec[-11] rec[-32]
DETECTOR(1, 2, 1) rec[-17] rec[-20] rec[-21] rec[-24] rec[-45]
DETECTOR(3, 2, 1) rec[-10] rec[-13] rec[-14] rec[-17] rec[-38]
DETECTOR(5, 2, 1) rec[-3] rec[-6] rec[-7] rec[-10] rec[-31]
DETECTOR(1, 4, 1) rec[-16] rec[-19] rec[-20] rec[-23] rec[-44]
DETECTOR(3, 4, 1) rec[-9] rec[-12] rec[-13] rec[-16] rec[-37]
DETECTOR(5, 4, 1) rec[-2] rec[-5] rec[-6] rec[-9] rec[-30]
DETECTOR(1, 6, 1) rec[-15] rec[-19] rec[-22] rec[-43]
DETECTOR(3, 6, 1) rec[-8] rec[-12] rec[-15] rec[-36]
DETECTOR(5, 6, 1) rec[-1] rec[-5] rec[-8] rec[-29]
OBSERVABLE_INCLUDE(0) rec[-22] rec[-23] rec[-24] rec[-25]
"""
In [ ]:
circuit = stim.Circuit(generate_circuit())
In [ ]:
circuit.without_noise().diagram('timeslice-svg')
In [ ]:
def count_logical_errors(circuit: stim.Circuit, num_shots: int) -> int:
    # Sample the circuit.
    sampler = circuit.compile_detector_sampler()
    sample =sampler.sample(num_shots, separate_observables=True)
    detection_events, observable_flips = sample
    detection_events = np.array(detection_events, order='C')

    # Configure a decoder using the circuit.
    detector_error_model = circuit.detector_error_model(decompose_errors=True)
    matcher = pymatching.Matching.from_detector_error_model(detector_error_model)

    # Run the decoder.
    predictions = []
    for i in range(num_shots):
        predictions.append(matcher.decode(detection_events[i]))
    predictions=np.array(predictions)

    # Count the mistakes.
    num_errors = 0
    for shot in range(num_shots):
        actual_for_shot = observable_flips[shot]
        predicted_for_shot = predictions[shot]
        if not np.array_equal(actual_for_shot, predicted_for_shot):
            num_errors += 1
    return num_errors
In [ ]:
circuit = stim.Circuit(generate_circuit(rounds=4))
num_shots = 100_000
num_logical_errors = count_logical_errors(circuit, num_shots)
print("there were", num_logical_errors, "wrong predictions (logical errors) out of", num_shots, "shots")
In [ ]:
dem = circuit.detector_error_model()
print(repr(dem))
In [ ]:
sampler = circuit.compile_sampler()
sample = sampler.sample(shots=1)

one_sample = sample[0]

for k in range(0, len(one_sample), 8):
    timeslice = one_sample[k:k+8]
    print("".join("1" if e else "_" for e in timeslice))
In [ ]:
circuit.without_noise().diagram(
    "detslice-with-ops-svg", 
    tick=range(0, 50),
)
In [ ]:
circuit.diagram("matchgraph-3d")
In [ ]:
surface_code_tasks = [
    sinter.Task(
        circuit = stim.Circuit(
            generate_circuit(
                rounds=3*4,
                after_clifford_depolarization=noise,
                after_reset_flip_probability=noise,
                before_measure_flip_probability=noise,
                before_round_data_depolarization=noise,
            )
        ),
        json_metadata={'d': d, 'r': d * 3, 'p': noise},
    )
    for d in [2]
    for noise in [0.0001,0.005,0.001,0.01,0.02,0.03,0.06,0.09,0.1,0.2,0.3]
]

collected_surface_code_stats: List[sinter.TaskStats] = sinter.collect(
    num_workers=4,
    tasks=surface_code_tasks,
    decoders=['pymatching'],
    max_shots=1_000_000,
    max_errors=5_000,
    print_progress=True,
)
In [ ]:
fig, ax = plt.subplots(1, 1)
sinter.plot_error_rate(
    ax=ax,
    stats=collected_surface_code_stats,
    x_func=lambda stat: stat.json_metadata['p'],
    group_func=lambda stat: stat.json_metadata['d'],
    failure_units_per_shot_func=lambda stat: stat.json_metadata['r'],
)
ax.loglog()
ax.set_title("Surface Code Error Rates per Round under Circuit Noise")
ax.set_xlabel("Phyical Error Rate")
ax.set_ylabel("Logical Error Rate per Round")
ax.grid(which='major')
ax.legend()
fig.set_dpi(120)  # Show it bigger
In [ ]:
from qccd import *

# Create architecture
arch = QCCDArch()

ionspacing = 5
trapspacing = 30
cooling_colour = 'red'
qubit_colour = 'lightblue'
traps = []

for i in range(3):
    ions = [QubitIon(qubit_colour, 'Q'), QubitIon(qubit_colour, 'Q'), QubitIon(qubit_colour, 'Q')]
    trap = arch.addManipulationTrap(x=i *trapspacing, y=0, ions=ions, color='grey', spacing=ionspacing, isHorizontal=True)
    traps.append(trap)

for t1, t2 in zip(traps[:-1], traps[1:]):
    arch.addEdge(t1, t2)

arch.refreshGraph()
fig, ax = plt.subplots()
# Display architecture
arch.display(fig,ax, showLabels=False)
In [ ]:
from qccd import *


# Create architecture
arch = QCCDArch()

ionspacing = 1
trapspacing = 6
qubit_colour = 'red'
qubit_colour = 'lightblue'
junction_colour = 'orange'
traps = []

ManipulationTrap.DEFAULT_SPACING = trapspacing
arch.SIZING = 0.5
ions11 = [QubitIon(qubit_colour, 'Q'), CoolingIon(qubit_colour, 'Q'), QubitIon(qubit_colour, 'Q'), CoolingIon(qubit_colour, 'Q')]
trap11 = arch.addManipulationTrap(x=0, y=0, ions=ions11, color='grey', spacing=ionspacing, isHorizontal=False)

ions12 = [QubitIon(qubit_colour, 'Q'), CoolingIon(qubit_colour, 'Q'), QubitIon(qubit_colour, 'Q'), CoolingIon(qubit_colour, 'Q')]
trap12 = arch.addManipulationTrap(x=trapspacing, y=0, ions=ions12, color='grey', spacing=ionspacing, isHorizontal=False)

junctionL = arch.addJunction(x=0, y=trapspacing, color=junction_colour)
junctionR = arch.addJunction(x=trapspacing, y=trapspacing, color=junction_colour)

ions21 = [QubitIon(qubit_colour, 'Q'), CoolingIon(qubit_colour, 'Q'), QubitIon(qubit_colour, 'Q'), CoolingIon(qubit_colour, 'Q')]
trap21 = arch.addManipulationTrap(x=0, y=2*trapspacing, ions=ions21, color='grey', spacing=ionspacing, isHorizontal=False)

ions22 = [QubitIon(qubit_colour, 'Q'), CoolingIon(qubit_colour, 'Q'), QubitIon(qubit_colour, 'Q'), CoolingIon(qubit_colour, 'Q')]
trap22 = arch.addManipulationTrap(x=trapspacing, y=2*trapspacing, ions=ions22, color='grey', spacing=ionspacing, isHorizontal=False)

crossing11 = arch.addEdge(trap11, junctionL)
arch.addEdge(trap12, junctionR)
crossing21 = arch.addEdge(trap21, junctionL)
arch.addEdge(trap22, junctionR)
arch.addEdge(junctionL, junctionR)

arch.refreshGraph()
# Display architecture


ops = (
    Split.physicalOperation(trap11, crossing11), 
    Move.physicalOperation(crossing11),
    JunctionCrossing.physicalOperation(junctionL, crossing11),
    JunctionCrossing.physicalOperation(junctionL, crossing21),
    Move.physicalOperation(crossing21),
    Merge.physicalOperation(trap21, crossing21),
    GateSwap.physicalOperation(trap=trap21, ion1=ions21[0], ion2=ions21[2]),
    TwoQubitMSGate.physicalOperation(ion1=ions22[0], ion2=ions22[2],trap=trap22),
    OneQubitGate.physicalOperation(ion=ions12[0], trap=trap12),
    Measurement.physicalOperation(ion=ions12[0], trap=trap12),
    QubitReset.physicalOperation(ion=ions12[0],trap=trap12))
parallelOps = paralleliseOperations(ops)

ncols = 2
fig, ax = plt.subplots()
arch.display(fig, ax, showLabels=False)
fig.set_size_inches(arch.WINDOW_SIZE[0]*0.5, arch.WINDOW_SIZE[1]*0.5)
fig,axs = plt.subplots(int(len(parallelOps)/2), 2)
axs = axs.flatten()
for i, (ax, parallelOp)in enumerate(zip(axs, parallelOps)):
    parallelOp.run()
    arch.refreshGraph()
    print(parallelOp.label+f' fidelity {parallelOp.fidelity()} operation time {parallelOp.operationTime()}')
    arch.display(fig=fig, ax=ax, title=parallelOp.label, operation=parallelOp, showLabels=False)

fig.set_size_inches(arch.WINDOW_SIZE[0]*0.5, arch.WINDOW_SIZE[1]*int(len(parallelOps)/2)*0.5)
In [ ]:
from qccd import *

noise = 1e-3
d=4
    
circuit = QCCDCircuit.generated(
    "surface_code:rotated_memory_z",
    rounds=1,
    distance=d,
    after_clifford_depolarization=noise,
    after_reset_flip_probability=noise,
    before_measure_flip_probability=noise,
    before_round_data_depolarization=noise,
)


rows = 4
cols = 4
capacity = 4

arch, instructions = circuit.processCircuit(rows=rows, cols=cols, trapCapacity=capacity, isHorizontal=(rows==1), clustering=1)

for i in instructions:
    print(i.label+str([ion.idx for ion in i.ions]))

capacities = [1,2]
rows = 6
cols = 6

label='part'
fig, axs = plt.subplots(1, len(capacities))
for ax, capacity in zip(axs, capacities):
    arch1, instructions = circuit.processCircuit(rows=rows, cols=cols, trapCapacity=capacity, isHorizontal=(rows==1), clustering=2)
    arch1.ION_SIZE = 2000

    for idx, (ion, pos) in circuit._ionMapping.items():
        ion.set(ion.idx, pos[0], pos[1])

    arch1.refreshGraph()
    edgesDups = []
    edgesPos = []
    ionsInvolved = set()
    score = 0
    for op in instructions:
        if not isinstance(op, TwoQubitMSGate):
            continue
        if not ionsInvolved.isdisjoint(op.ions):
            score+=1
            ionsInvolved=set()
        edgesDups.append((op.ions, score))
        edgesPos.append((op.ionsActedIdxs, score))
        ionsInvolved=ionsInvolved.union(op.ions)
    scores = [score for ((ion1,ion2), score) in edgesDups]
    color_from_score = {s: (5+i*i)*2 for i, s in enumerate( sorted(list(set(scores)), reverse=True))}
    arch1._manipulationTraps.append(([(ion1.idx, ion2.idx) for ((ion1,ion2), score) in edgesDups], [color_from_score[score] for ((ion1,ion2), score) in edgesDups]))
    arch1.display(fig, ax, showLabels=False, showEdges=False, show_junction=False)
    ax.set_title(f"{label} {capacities}")


circuit.without_noise().diagram('timeline-svg')
QUBIT_RESET[2]
QUBIT_RESET[3]
QUBIT_RESET[15]
QUBIT_RESET[16]
QUBIT_RESET[26]
QUBIT_RESET[27]
QUBIT_RESET[39]
QUBIT_RESET[6]
QUBIT_RESET[19]
QUBIT_RESET[30]
QUBIT_RESET[42]
QUBIT_RESET[9]
QUBIT_RESET[33]
QUBIT_RESET[34]
QUBIT_RESET[45]
QUBIT_RESET[12]
QUBIT_RESET[1]
QUBIT_RESET[14]
QUBIT_RESET[25]
QUBIT_RESET[38]
QUBIT_RESET[5]
QUBIT_RESET[18]
QUBIT_RESET[29]
QUBIT_RESET[41]
QUBIT_RESET[8]
QUBIT_RESET[21]
QUBIT_RESET[32]
QUBIT_RESET[44]
QUBIT_RESET[11]
QUBIT_RESET[23]
QUBIT_RESET[36]
RY[1]
RX[1]
RY[14]
RX[14]
RY[38]
RX[38]
RY[29]
RX[29]
RY[8]
RX[8]
RY[44]
RX[44]
RY[23]
RX[23]
RY[36]
RX[36]
RY[1]
RX[1]
RX[3]
TWO_QUBIT_MS_GATE[1, 3]
RY[1]
RY[29]
RX[29]
RX[30]
TWO_QUBIT_MS_GATE[29, 30]
RY[29]
RY[38]
RX[38]
RX[39]
TWO_QUBIT_MS_GATE[38, 39]
RY[38]
RY[44]
RX[44]
RX[45]
TWO_QUBIT_MS_GATE[44, 45]
RY[44]
RY[14]
RX[14]
RX[16]
TWO_QUBIT_MS_GATE[14, 16]
RY[14]
RY[8]
RX[8]
RX[9]
TWO_QUBIT_MS_GATE[8, 9]
RY[8]
RY[19]
RX[19]
RX[18]
TWO_QUBIT_MS_GATE[19, 18]
RY[19]
RY[27]
RX[27]
RX[25]
TWO_QUBIT_MS_GATE[27, 25]
RY[27]
RY[34]
RX[34]
RX[32]
TWO_QUBIT_MS_GATE[34, 32]
RY[34]
RY[42]
RX[42]
RX[41]
TWO_QUBIT_MS_GATE[42, 41]
RY[42]
RY[6]
RX[6]
RX[5]
TWO_QUBIT_MS_GATE[6, 5]
RY[6]
RY[12]
RX[12]
RX[11]
TWO_QUBIT_MS_GATE[12, 11]
RY[12]
RY[1]
RX[1]
RX[2]
TWO_QUBIT_MS_GATE[1, 2]
RY[1]
RY[29]
RX[29]
RX[19]
TWO_QUBIT_MS_GATE[29, 19]
RY[29]
RY[38]
RX[38]
RX[27]
TWO_QUBIT_MS_GATE[38, 27]
RY[38]
RY[44]
RX[44]
RX[34]
TWO_QUBIT_MS_GATE[44, 34]
RY[44]
RY[14]
RX[14]
RX[15]
TWO_QUBIT_MS_GATE[14, 15]
RY[14]
RY[8]
RX[8]
RX[42]
TWO_QUBIT_MS_GATE[8, 42]
RY[8]
RY[26]
RX[26]
RX[18]
TWO_QUBIT_MS_GATE[26, 18]
RY[26]
RY[3]
RX[3]
RX[25]
TWO_QUBIT_MS_GATE[3, 25]
RY[3]
RY[30]
RX[30]
RX[32]
TWO_QUBIT_MS_GATE[30, 32]
RY[30]
RY[39]
RX[39]
RX[41]
TWO_QUBIT_MS_GATE[39, 41]
RY[39]
RY[16]
RX[16]
RX[5]
TWO_QUBIT_MS_GATE[16, 5]
RY[16]
RY[9]
RX[9]
RX[11]
TWO_QUBIT_MS_GATE[9, 11]
RY[9]
RY[29]
RX[29]
RX[27]
TWO_QUBIT_MS_GATE[29, 27]
RY[29]
RY[23]
RX[23]
RX[34]
TWO_QUBIT_MS_GATE[23, 34]
RY[23]
RY[38]
RX[38]
RX[15]
TWO_QUBIT_MS_GATE[38, 15]
RY[38]
RY[44]
RX[44]
RX[42]
TWO_QUBIT_MS_GATE[44, 42]
RY[44]
RY[8]
RX[8]
RX[6]
TWO_QUBIT_MS_GATE[8, 6]
RY[8]
RY[36]
RX[36]
RX[12]
TWO_QUBIT_MS_GATE[36, 12]
RY[36]
RY[26]
RX[26]
RX[25]
TWO_QUBIT_MS_GATE[26, 25]
RY[26]
RY[33]
RX[33]
RX[32]
TWO_QUBIT_MS_GATE[33, 32]
RY[33]
RY[30]
RX[30]
RX[41]
TWO_QUBIT_MS_GATE[30, 41]
RY[30]
RY[39]
RX[39]
RX[5]
TWO_QUBIT_MS_GATE[39, 5]
RY[39]
RY[45]
RX[45]
RX[11]
TWO_QUBIT_MS_GATE[45, 11]
RY[45]
RY[9]
RX[9]
RX[21]
TWO_QUBIT_MS_GATE[9, 21]
RY[9]
RY[29]
RX[29]
RX[26]
TWO_QUBIT_MS_GATE[29, 26]
RY[29]
RY[23]
RX[23]
RX[33]
TWO_QUBIT_MS_GATE[23, 33]
RY[23]
RY[38]
RX[38]
RX[3]
TWO_QUBIT_MS_GATE[38, 3]
RY[38]
RY[44]
RX[44]
RX[30]
TWO_QUBIT_MS_GATE[44, 30]
RY[44]
RY[8]
RX[8]
RX[39]
TWO_QUBIT_MS_GATE[8, 39]
RY[8]
RY[36]
RX[36]
RX[45]
TWO_QUBIT_MS_GATE[36, 45]
RY[36]
RY[2]
RX[2]
RX[25]
TWO_QUBIT_MS_GATE[2, 25]
RY[2]
RY[19]
RX[19]
RX[32]
TWO_QUBIT_MS_GATE[19, 32]
RY[19]
RY[27]
RX[27]
RX[41]
TWO_QUBIT_MS_GATE[27, 41]
RY[27]
RY[15]
RX[15]
RX[5]
TWO_QUBIT_MS_GATE[15, 5]
RY[15]
RY[42]
RX[42]
RX[11]
TWO_QUBIT_MS_GATE[42, 11]
RY[42]
RY[6]
RX[6]
RX[21]
TWO_QUBIT_MS_GATE[6, 21]
RY[6]
RY[1]
RX[1]
RY[14]
RX[14]
RY[38]
RX[38]
RY[29]
RX[29]
RY[8]
RX[8]
RY[44]
RX[44]
RY[23]
RX[23]
RY[36]
RX[36]
MEASUREMENT[1]
MEASUREMENT[14]
MEASUREMENT[25]
MEASUREMENT[38]
MEASUREMENT[5]
MEASUREMENT[18]
MEASUREMENT[29]
MEASUREMENT[41]
MEASUREMENT[8]
MEASUREMENT[21]
MEASUREMENT[32]
MEASUREMENT[44]
MEASUREMENT[11]
MEASUREMENT[23]
MEASUREMENT[36]
QUBIT_RESET[1]
QUBIT_RESET[14]
QUBIT_RESET[25]
QUBIT_RESET[38]
QUBIT_RESET[5]
QUBIT_RESET[18]
QUBIT_RESET[29]
QUBIT_RESET[41]
QUBIT_RESET[8]
QUBIT_RESET[21]
QUBIT_RESET[32]
QUBIT_RESET[44]
QUBIT_RESET[11]
QUBIT_RESET[23]
QUBIT_RESET[36]
MEASUREMENT[2]
MEASUREMENT[3]
MEASUREMENT[15]
MEASUREMENT[16]
MEASUREMENT[26]
MEASUREMENT[27]
MEASUREMENT[39]
MEASUREMENT[6]
MEASUREMENT[19]
MEASUREMENT[30]
MEASUREMENT[42]
MEASUREMENT[9]
MEASUREMENT[33]
MEASUREMENT[34]
MEASUREMENT[45]
MEASUREMENT[12]
Out[ ]:
No description has been provided for this image
No description has been provided for this image
In [ ]:
from qccd import *

noise = 1e-3
d=4
    
circuit = QCCDCircuit.generated(
    "surface_code:rotated_memory_z",
    rounds=1,
    distance=d,
    after_clifford_depolarization=noise,
    after_reset_flip_probability=noise,
    before_measure_flip_probability=noise,
    before_round_data_depolarization=noise,
)


rows = 4
cols = 4
capacity = 4

arch, instructions = circuit.processCircuit(rows=rows, cols=cols, trapCapacity=capacity, isHorizontal=(rows==1), clustering=1)

for i in instructions:
    print(i.label+str([ion.idx for ion in i.ions]))

capacities = [2,3,4,5,6,7]
rows = 6
cols = 6


for i, label in enumerate(('basic', 'pos', 'part')):
    fig, axs = plt.subplots(1, len(capacities))
    for ax, capacity in zip(axs, capacities):
        arch1, instructions = circuit.processCircuit(rows=rows, cols=cols, trapCapacity=capacity, isHorizontal=(rows==1), clustering=i)
        arch1.ION_SIZE = 2000

        for idx, (ion, pos) in circuit._ionMapping.items():
            ion.set(ion.idx, pos[0], pos[1])

        arch1.refreshGraph()
        edgesDups = []
        edgesPos = []
        ionsInvolved = set()
        score = 0
        for op in instructions:
            if not isinstance(op, TwoQubitMSGate):
                continue
            if not ionsInvolved.isdisjoint(op.ions):
                score+=1
                ionsInvolved=set()
            edgesDups.append((op.ions, score))
            edgesPos.append((op.ionsActedIdxs, score))
            ionsInvolved=ionsInvolved.union(op.ions)
        scores = [score for ((ion1,ion2), score) in edgesDups]
        color_from_score = {s: (5+i*i)*0 for i, s in enumerate( sorted(list(set(scores)), reverse=True))}
        arch1._manipulationTraps.append(([(ion1.idx, ion2.idx) for ((ion1,ion2), score) in edgesDups], [color_from_score[score] for ((ion1,ion2), score) in edgesDups]))
        arch1.display(fig, ax, showLabels=False, showEdges=False, show_junction=False)
        ax.set_title(f"{label} {capacities}")


circuit.without_noise().diagram('timeline-svg')
QUBIT_RESET[2]
QUBIT_RESET[3]
QUBIT_RESET[15]
QUBIT_RESET[16]
QUBIT_RESET[26]
QUBIT_RESET[27]
QUBIT_RESET[39]
QUBIT_RESET[6]
QUBIT_RESET[19]
QUBIT_RESET[30]
QUBIT_RESET[42]
QUBIT_RESET[9]
QUBIT_RESET[33]
QUBIT_RESET[34]
QUBIT_RESET[45]
QUBIT_RESET[12]
QUBIT_RESET[1]
QUBIT_RESET[14]
QUBIT_RESET[25]
QUBIT_RESET[38]
QUBIT_RESET[5]
QUBIT_RESET[18]
QUBIT_RESET[29]
QUBIT_RESET[41]
QUBIT_RESET[8]
QUBIT_RESET[21]
QUBIT_RESET[32]
QUBIT_RESET[44]
QUBIT_RESET[11]
QUBIT_RESET[23]
QUBIT_RESET[36]
RY[1]
RX[1]
RY[14]
RX[14]
RY[38]
RX[38]
RY[29]
RX[29]
RY[8]
RX[8]
RY[44]
RX[44]
RY[23]
RX[23]
RY[36]
RX[36]
RY[1]
RX[1]
RX[3]
TWO_QUBIT_MS_GATE[1, 3]
RY[1]
RY[29]
RX[29]
RX[30]
TWO_QUBIT_MS_GATE[29, 30]
RY[29]
RY[38]
RX[38]
RX[39]
TWO_QUBIT_MS_GATE[38, 39]
RY[38]
RY[44]
RX[44]
RX[45]
TWO_QUBIT_MS_GATE[44, 45]
RY[44]
RY[14]
RX[14]
RX[16]
TWO_QUBIT_MS_GATE[14, 16]
RY[14]
RY[8]
RX[8]
RX[9]
TWO_QUBIT_MS_GATE[8, 9]
RY[8]
RY[19]
RX[19]
RX[18]
TWO_QUBIT_MS_GATE[19, 18]
RY[19]
RY[27]
RX[27]
RX[25]
TWO_QUBIT_MS_GATE[27, 25]
RY[27]
RY[34]
RX[34]
RX[32]
TWO_QUBIT_MS_GATE[34, 32]
RY[34]
RY[42]
RX[42]
RX[41]
TWO_QUBIT_MS_GATE[42, 41]
RY[42]
RY[6]
RX[6]
RX[5]
TWO_QUBIT_MS_GATE[6, 5]
RY[6]
RY[12]
RX[12]
RX[11]
TWO_QUBIT_MS_GATE[12, 11]
RY[12]
RY[1]
RX[1]
RX[2]
TWO_QUBIT_MS_GATE[1, 2]
RY[1]
RY[29]
RX[29]
RX[19]
TWO_QUBIT_MS_GATE[29, 19]
RY[29]
RY[38]
RX[38]
RX[27]
TWO_QUBIT_MS_GATE[38, 27]
RY[38]
RY[44]
RX[44]
RX[34]
TWO_QUBIT_MS_GATE[44, 34]
RY[44]
RY[14]
RX[14]
RX[15]
TWO_QUBIT_MS_GATE[14, 15]
RY[14]
RY[8]
RX[8]
RX[42]
TWO_QUBIT_MS_GATE[8, 42]
RY[8]
RY[26]
RX[26]
RX[18]
TWO_QUBIT_MS_GATE[26, 18]
RY[26]
RY[3]
RX[3]
RX[25]
TWO_QUBIT_MS_GATE[3, 25]
RY[3]
RY[30]
RX[30]
RX[32]
TWO_QUBIT_MS_GATE[30, 32]
RY[30]
RY[39]
RX[39]
RX[41]
TWO_QUBIT_MS_GATE[39, 41]
RY[39]
RY[16]
RX[16]
RX[5]
TWO_QUBIT_MS_GATE[16, 5]
RY[16]
RY[9]
RX[9]
RX[11]
TWO_QUBIT_MS_GATE[9, 11]
RY[9]
RY[29]
RX[29]
RX[27]
TWO_QUBIT_MS_GATE[29, 27]
RY[29]
RY[23]
RX[23]
RX[34]
TWO_QUBIT_MS_GATE[23, 34]
RY[23]
RY[38]
RX[38]
RX[15]
TWO_QUBIT_MS_GATE[38, 15]
RY[38]
RY[44]
RX[44]
RX[42]
TWO_QUBIT_MS_GATE[44, 42]
RY[44]
RY[8]
RX[8]
RX[6]
TWO_QUBIT_MS_GATE[8, 6]
RY[8]
RY[36]
RX[36]
RX[12]
TWO_QUBIT_MS_GATE[36, 12]
RY[36]
RY[26]
RX[26]
RX[25]
TWO_QUBIT_MS_GATE[26, 25]
RY[26]
RY[33]
RX[33]
RX[32]
TWO_QUBIT_MS_GATE[33, 32]
RY[33]
RY[30]
RX[30]
RX[41]
TWO_QUBIT_MS_GATE[30, 41]
RY[30]
RY[39]
RX[39]
RX[5]
TWO_QUBIT_MS_GATE[39, 5]
RY[39]
RY[45]
RX[45]
RX[11]
TWO_QUBIT_MS_GATE[45, 11]
RY[45]
RY[9]
RX[9]
RX[21]
TWO_QUBIT_MS_GATE[9, 21]
RY[9]
RY[29]
RX[29]
RX[26]
TWO_QUBIT_MS_GATE[29, 26]
RY[29]
RY[23]
RX[23]
RX[33]
TWO_QUBIT_MS_GATE[23, 33]
RY[23]
RY[38]
RX[38]
RX[3]
TWO_QUBIT_MS_GATE[38, 3]
RY[38]
RY[44]
RX[44]
RX[30]
TWO_QUBIT_MS_GATE[44, 30]
RY[44]
RY[8]
RX[8]
RX[39]
TWO_QUBIT_MS_GATE[8, 39]
RY[8]
RY[36]
RX[36]
RX[45]
TWO_QUBIT_MS_GATE[36, 45]
RY[36]
RY[2]
RX[2]
RX[25]
TWO_QUBIT_MS_GATE[2, 25]
RY[2]
RY[19]
RX[19]
RX[32]
TWO_QUBIT_MS_GATE[19, 32]
RY[19]
RY[27]
RX[27]
RX[41]
TWO_QUBIT_MS_GATE[27, 41]
RY[27]
RY[15]
RX[15]
RX[5]
TWO_QUBIT_MS_GATE[15, 5]
RY[15]
RY[42]
RX[42]
RX[11]
TWO_QUBIT_MS_GATE[42, 11]
RY[42]
RY[6]
RX[6]
RX[21]
TWO_QUBIT_MS_GATE[6, 21]
RY[6]
RY[1]
RX[1]
RY[14]
RX[14]
RY[38]
RX[38]
RY[29]
RX[29]
RY[8]
RX[8]
RY[44]
RX[44]
RY[23]
RX[23]
RY[36]
RX[36]
MEASUREMENT[1]
MEASUREMENT[14]
MEASUREMENT[25]
MEASUREMENT[38]
MEASUREMENT[5]
MEASUREMENT[18]
MEASUREMENT[29]
MEASUREMENT[41]
MEASUREMENT[8]
MEASUREMENT[21]
MEASUREMENT[32]
MEASUREMENT[44]
MEASUREMENT[11]
MEASUREMENT[23]
MEASUREMENT[36]
QUBIT_RESET[1]
QUBIT_RESET[14]
QUBIT_RESET[25]
QUBIT_RESET[38]
QUBIT_RESET[5]
QUBIT_RESET[18]
QUBIT_RESET[29]
QUBIT_RESET[41]
QUBIT_RESET[8]
QUBIT_RESET[21]
QUBIT_RESET[32]
QUBIT_RESET[44]
QUBIT_RESET[11]
QUBIT_RESET[23]
QUBIT_RESET[36]
MEASUREMENT[2]
MEASUREMENT[3]
MEASUREMENT[15]
MEASUREMENT[16]
MEASUREMENT[26]
MEASUREMENT[27]
MEASUREMENT[39]
MEASUREMENT[6]
MEASUREMENT[19]
MEASUREMENT[30]
MEASUREMENT[42]
MEASUREMENT[9]
MEASUREMENT[33]
MEASUREMENT[34]
MEASUREMENT[45]
MEASUREMENT[12]
Out[ ]:
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [ ]:
from qccd_circuit import *

noise=0.01
circuit = QCCDCircuit.generated(
    "surface_code:rotated_memory_z",
    rounds=1,
    distance=4,
    after_clifford_depolarization=noise,
    after_reset_flip_probability=noise,
    before_measure_flip_probability=noise,
    before_round_data_depolarization=noise,
)
circuit.without_noise().diagram('timeslice-svg')
Out[ ]:
No description has been provided for this image
In [ ]:
arch, instructions = circuit.processCircuitAugmentedGrid(rows=3, cols=3, trapCapacity=7)

oldPositions = {}
for idx, (ion, pos) in circuit._ionMapping.items():
    oldPositions[idx] = ion.pos
    ion.set(ion.idx, pos[0], pos[1], parent=ion.parent)

fig,ax =plt.subplots()
arch1 = arch
arch1.refreshGraph()
edgesDups = []
edgesPos = []
ionsInvolved = set()
score = 0
for op in instructions:
    if not isinstance(op, TwoQubitMSGate):
        continue
    if not ionsInvolved.isdisjoint(op.ions):
        score+=1
        ionsInvolved=set()
    edgesDups.append((op.ions, score))
    edgesPos.append((op.ionsActedIdxs, score))
    ionsInvolved=ionsInvolved.union(op.ions)
scores = [score for ((ion1,ion2), score) in edgesDups]
color_from_score = {s: (5+i*i)*0 for i, s in enumerate( sorted(list(set(scores)), reverse=True))}
arch1._manipulationTraps.append(([(ion1.idx, ion2.idx) for ((ion1,ion2), score) in edgesDups], [color_from_score[score] for ((ion1,ion2), score) in edgesDups]))
arch1.display(fig, ax, showLabels=False, showEdges=False, show_junction=False)
fig.set_size_inches(arch.WINDOW_SIZE[0]*0.6, arch.WINDOW_SIZE[1]*0.6)
arch1._manipulationTraps = arch1._manipulationTraps[:-1]
for idx, (ion, pos) in circuit._ionMapping.items():
    ion.set(ion.idx, oldPositions[idx][0], oldPositions[idx][1], parent=ion.parent)


arch.refreshGraph()

fig,ax =plt.subplots()
arch.display(fig, ax, showLabels=False)
fig.set_size_inches(arch.WINDOW_SIZE[0]*0.8, arch.WINDOW_SIZE[1]*0.8)

allOps, _ = arch.processOperationsViaForwarding(instructions)
parallelOpsMap = paralleliseOperations(allOps)
parallelOps = list(dict(parallelOpsMap).values())

# arch = circuit.resetArch()
# arch.refreshGraph()


# axsPerFig = 2
# figs = []
# axs = []
# numFigs = int(len(parallelOps)/(axsPerFig))+int((len(parallelOps)%(axsPerFig))>0)
# for _ in range(numFigs):
#     fig, axs_ = plt.subplots(1, axsPerFig, sharex=False, sharey=False)
#     figs.extend([fig]*axsPerFig)
#     axs.extend(axs_)

# figs = figs[:len(parallelOps)]
# axs = axs[:len(parallelOps)]

# for fig, ax, parallelOp in zip(figs, axs, parallelOps):
#     arch.display(fig, ax, operation=parallelOp, runOps=True, showLabels=False)
  
print(f"total number of qubit operations: {len(instructions)}")
print(f"total number of operations: {len(allOps)}")
print(f"time for operations: {max(parallelOpsMap.keys())}")
total number of qubit operations: 349
total number of operations: 610
time for operations: 0.02321800000000003
No description has been provided for this image
No description has been provided for this image
In [ ]:
from qccd_circuit import *

noise=0.01
circuit = QCCDCircuit.generated(
    "surface_code:rotated_memory_z",
    rounds=1,
    distance=4,
    after_clifford_depolarization=noise,
    after_reset_flip_probability=noise,
    before_measure_flip_probability=noise,
    before_round_data_depolarization=noise,
)
arch, instructions = circuit.processCircuitAugmentedGrid(rows=5, cols=5, trapCapacity=1)

oldPositions = {}
for idx, (ion, pos) in circuit._ionMapping.items():
    oldPositions[idx] = ion.pos
    ion.set(ion.idx, pos[0], pos[1], parent=ion.parent)

fig,ax =plt.subplots()
arch1 = arch
arch1.refreshGraph()
edgesDups = []
edgesPos = []
ionsInvolved = set()
score = 0
for op in instructions:
    if not isinstance(op, TwoQubitMSGate):
        continue
    if not ionsInvolved.isdisjoint(op.ions):
        score+=1
        ionsInvolved=set()
    edgesDups.append((op.ions, score))
    edgesPos.append((op.ionsActedIdxs, score))
    ionsInvolved=ionsInvolved.union(op.ions)
scores = [score for ((ion1,ion2), score) in edgesDups]
color_from_score = {s: (5+i*i)*0 for i, s in enumerate( sorted(list(set(scores)), reverse=True))}
arch1._manipulationTraps.append(([(ion1.idx, ion2.idx) for ((ion1,ion2), score) in edgesDups], [color_from_score[score] for ((ion1,ion2), score) in edgesDups]))
arch1.display(fig, ax, showLabels=False, showEdges=False, show_junction=False)
fig.set_size_inches(arch.WINDOW_SIZE[0]*0.5, arch.WINDOW_SIZE[1]*0.5)
arch1._manipulationTraps = arch1._manipulationTraps[:-1]
for idx, (ion, pos) in circuit._ionMapping.items():
    ion.set(ion.idx, oldPositions[idx][0], oldPositions[idx][1], parent=ion.parent)

arch.refreshGraph()

fig,ax =plt.subplots()
arch.display(fig, ax, title='map complete', showLabels=False)
fig.set_size_inches(arch.WINDOW_SIZE[0]*0.5, arch.WINDOW_SIZE[1]*0.5)

allOps, barriers = arch.processOperationsViaForwarding(instructions)
parallelOpsMap = paralleliseOperationsWithBarriers(allOps, barriers)
parallelOps = list(dict(parallelOpsMap).values())

circuit.simulate(allOps)

arch = circuit.resetArch()
arch.refreshGraph()


print(f"total number of qubit operations: {len(instructions)}")
print(f"total number of operations: {len(allOps)}")
print(f"time for operations: {max(parallelOpsMap.keys())}")
total number of qubit operations: 349
total number of operations: 955
time for operations: 0.012392
No description has been provided for this image
No description has been provided for this image
In [ ]:
print(len(paralleliseOperationsSimple(allOps)))
In [ ]:
parallelOps = paralleliseOperationsSimple(allOps)
arch = circuit.resetArch()
arch.refreshGraph()


axsPerFig = 2
figs = []
axs = []
numFigs = int(len(parallelOps)/(axsPerFig))+int((len(parallelOps)%(axsPerFig))>0)
for _ in range(numFigs):
    fig, axs_ = plt.subplots(1, axsPerFig, sharex=False, sharey=False)
    figs.extend([fig]*axsPerFig)
    axs.extend(axs_)

figs = figs[:len(parallelOps)]
axs = axs[:len(parallelOps)]

for fig, ax, parallelOp in zip(figs, axs, parallelOps):
    arch.display(fig, ax, operation=parallelOp, runOps=True, showLabels=False)
  
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
%matplotlib inline
from IPython.display import HTML, clear_output, display

arch = circuit.resetArch()
arch.refreshGraph()
fig,ax=plt.subplots()

parallelOpsSimple = paralleliseOperationsSimple(allOps)
def update(frame, _arch):
    ax.clear()
    # Clear the output for dynamic display in Jupyter notebook
    clear_output(wait=True)
    if frame>0:
        op: ParallelOperation = parallelOpsSimple[frame-1]
        title = f"Operation: {op.label}"
        _arch.display(fig, ax, title, operation=op, runOps=True, showLabels=False)
    else:
        _arch.display(fig, ax, showLabels=False)
    return display(fig)
import time
time.sleep(10)
ani = FuncAnimation(fig, lambda frame: update(frame, arch), frames=len(parallelOpsSimple)+1, repeat=False)
ani.save('animation.mp4', writer='ffmpeg', fps=10, dpi=100)
# HTML(ani.to_jshtml())
No description has been provided for this image
No description has been provided for this image
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
%matplotlib inline
from IPython.display import HTML, clear_output, display

arch = circuit.resetArch()
arch.refreshGraph()
fig,ax=plt.subplots()
parallelOpsTimes = sorted(parallelOpsMap.keys())


def update(frame, _arch):
    ax.clear()
    # Clear the output for dynamic display in Jupyter notebook
    clear_output(wait=True)
    if frame>0:
        op: ParallelOperation = parallelOps[frame-1]
        title = f"Operation {round(parallelOpsTimes[frame-1],6)}: {op.label}"
        _arch.display(fig, ax, title, operation=op, runOps=True, showLabels=False)
    else:
        _arch.display(fig, ax, showLabels=False)
    return display(fig)
import time
time.sleep(10)
ani = FuncAnimation(fig, lambda frame: update(frame, arch), frames=len(parallelOps)+1, repeat=False)
ani.save('animation.mp4', writer='ffmpeg', fps=10, dpi=100)
HTML(ani.to_jshtml())
No description has been provided for this image
No description has been provided for this image
In [ ]:
# gate improvement that is expected in the next few years is about 10X or maximum 100X, not more than that
gate_improvements = [0.025, 0.05, 0.1, 0.25, 0.5, 1.0, 5.0, 10.0, 12.5, 15.0, 20.0, 25.0, 35.0, 50.0, 100.0]
distances = [3,5,7,9,11]
# capacities refer to maximum number of qubits we place in a single trap for the initial QCCD configuration
# true maximum trap capacity, i.e. the maximum number of ions each trap can hold in any instance, is 5 times this initial configuration maximum value 
capacities = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,30,60,120]


data = {k: {"Forwarding": {d: {c: None for c in capacities} for d in distances}, "Routing": {d: {c: None for c in capacities} for d in distances}} for k in ("PhysicalErrorRates", "Operations", "ElapsedTime", "MeanConcurrency", "QubitOperations", "LogicalErrorRates")}
In [ ]:
for k1 in data.keys():
    for k2 in data[k1].keys():
        for k3 in data[k1][k2].keys():
            print(k1, k2, k3, data[k1][k2][k3])
In [ ]:
data = ...
In [ ]:
with open('saved_Data.txt', 'w') as f:
    f.write(str(data))
In [ ]:
from qccd_circuit import *

num_shots = 1_000_000
data: Dict[str, Dict[str, Dict[int, Dict[int, Any]]]]


import concurrent.futures
import logging

import os

num_cores = os.cpu_count()
num_processes = num_cores 

logger = logging.getLogger()
logger.setLevel(logging.INFO)
handler = logging.FileHandler("process_log.txt")
formatter = logging.Formatter('%(processName)s - %(levelname)s - %(message)s')
handler.setFormatter(formatter)
logger.addHandler(handler)

with concurrent.futures.ProcessPoolExecutor(max_workers=num_processes) as executor:
    futures = []
    for d in distances:
        for c in capacities:
            logger.info(f"Submitting task for distance {d} and capacity {c}.")
            futures.append(executor.submit(process_circuit, d, c, gate_improvements, num_shots))

    for future in concurrent.futures.as_completed(futures):
        if future.exception() is not None:
            logger.error(f"Exception Occured for future", exc_info=future.exception())
            continue

        try:
            result = future.result()

            d = result["Distance"]
            c = result["Capacity"]
            for label in result["ElapsedTime"]:
                data["ElapsedTime"][label][d][c] = result["ElapsedTime"][label]
                data["Operations"][label][d][c] = result["Operations"][label]
                data["MeanConcurrency"][label][d][c] = result["MeanConcurrency"][label]
                data["QubitOperations"][label][d][c] = result["QubitOperations"][label]
                data["LogicalErrorRates"][label][d][c] = result["LogicalErrorRates"][label]
                data["PhysicalErrorRates"][label][d][c] = result["PhysicalErrorRates"][label]

            
            logger.info(f"Processing result for distance {result['Distance']}, capacity {result['Capacity']}.")
        except Exception as e:
            logger.error(f"Exception Occured for future", exc_info=True)
            continue

        
In [ ]:
for k1 in data.keys():
    for k2 in data[k1].keys():
        for k3 in data[k1][k2].keys():
            data[k1][k2][k3] = [data[k1][k2][k3][c] for c in capacities]


data: Dict[str, Dict[str, Dict[int, Dict[int, Any]]]]
In [ ]:
import matplotlib as mpl
import matplotlib.colors as mplcolors
In [ ]:
cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, 3))
colors = [mplcolors.to_hex(c) for c in colors]

fig1, ax1 = plt.subplots()

iLabel, label = 0, "Forwarding"
ax1.plot(distances, [data["Operations"][label][d][0] for d in distances], c=colors[0], linestyle=("dashed" if iLabel==1 else "solid"), label=f"Total Operations for QEC shot")
ax1.plot(distances, [int(data["Operations"][label][d][0]/data["MeanConcurrency"][label][d][0]) for d in distances], c=colors[1], linestyle=("dashed" if iLabel==1 else "solid"), label=f"Mean Circuit Depth for QEC shot")
ax1.scatter(distances, [data["Operations"][label][d][0] for d in distances], c=colors[0], marker=("v" if iLabel==1 else "o"))
ax1.scatter(distances, [int(data["Operations"][label][d][0]/data["MeanConcurrency"][label][d][0]) for d in distances], c=colors[1], marker=("v" if iLabel==1 else "o"))
ax1.legend()
ax1.set_title('Importance of Effective Scheduling')
ax1.set_xticks(distances)
ax1.set_xlabel('Distances')
ax1.set_ylabel('Num. Ops.')
ax1.grid(which='major')
In [ ]:
_distances = [3,5,7,9,11]

cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, 3))
colors = [mplcolors.to_hex(c) for c in colors]

fig1, ax1 = plt.subplots()

iLabel, label = 0, "Forwarding"
ax1.plot(_distances, [data["Operations"][label][d][0] for d in _distances], c=colors[0], linestyle=("dashed" if iLabel==1 else "solid"), label=f"Total Number of Operations")
ax1.plot(_distances, [data["QubitOperations"][label][d][0] for d in _distances], c=colors[1], linestyle=("dashed" if iLabel==1 else "solid"), label=f"Number of Gates")
ax1.scatter(_distances, [data["Operations"][label][d][0] for d in _distances], c=colors[0], marker=("v" if iLabel==1 else "o"))
ax1.scatter(_distances, [data["QubitOperations"][label][d][0] for d in _distances], c=colors[1], marker=("v" if iLabel==1 else "o"))
ax1.legend()
ax1.set_title('Importance of Effective Mapping')
ax1.set_xlabel('Distances')
ax1.set_xticks(_distances)
ax1.set_ylabel('Operations')
ax1.grid(which='major')
In [ ]:
cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, len(distances)+1))
colors = [mplcolors.to_hex(c) for c in colors]

fig1, ax1 = plt.subplots()

iLabel, label = 0, "Forwarding"
for j, d in enumerate(distances):
    ax1.plot(capacities, [data["Operations"][label][d][i] for i in range(len(capacities))], c=colors[j], linestyle=("dashed" if iLabel==1 else "solid"), label=f"d={d}")
    ax1.scatter(capacities, [data["Operations"][label][d][i] for i in range(len(capacities))], c=colors[j], marker=("v" if iLabel==1 else "o"))

ax1.legend()
ax1.set_title('Total Number of Operations in One Shot of a QEC Cycle')
ax1.set_xlabel('Trap Capacity')
ax1.set_ylabel('Num. Ops.')
ax1.grid(which='major')



fig1, ax1 = plt.subplots()

iLabel, label = 0, "Forwarding"
for j, d in enumerate([distances[0], distances[-1]]):
    ax1.plot(capacities, [data["MeanConcurrency"][label][d][i] for i in range(len(capacities))], c=colors[j], linestyle=("dashed" if iLabel==1 else "solid"), label=f"d={d}")
    ax1.scatter(capacities, [data["MeanConcurrency"][label][d][i] for i in range(len(capacities))], c=colors[j], marker=("v" if iLabel==1 else "o"))

ax1.legend()
ax1.set_title('Concurrent Operations in a Timeslice of a QEC Cycle')
ax1.set_xlabel('Trap Capacity')
ax1.set_ylabel('Num. Ops.')
ax1.grid(which='major')
In [ ]:
cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, len(capacities)+1))
colors = [mplcolors.to_hex(c) for c in colors]

fig5, ax5 = plt.subplots()
for i, c in enumerate(capacities):
    if i%4==0:
        ax5.plot(distances, [round(data["MeanConcurrency"][label][d][i]/(d**2-1), 3) for d in distances], c=colors[i], linestyle=("dashed" if iLabel==1 else "solid"), label=f"c={c}")
        ax5.scatter(distances, [round(data["MeanConcurrency"][label][d][i]/(d**2-1), 3) for d in distances], c=colors[i], marker=("v" if iLabel==1 else "o"))

ax5.legend()
ax5.set_title('Concurrency vs Capacity')
ax5.set_xticks(distances)
ax5.set_xlabel('Code Distance')
ax5.set_ylabel('Concurrent Operations per Plaquette')
ax5.grid(which='major')


fig4, ax4 = plt.subplots()
for i, c in enumerate(capacities):
    if i%4==0:
        ax4.plot(distances, [data["ElapsedTime"][label][d][i] for d in distances], c=colors[i], linestyle=("dashed" if iLabel==1 else "solid"),label=f"c={c}")
        ax4.scatter(distances, [data["ElapsedTime"][label][d][i] for d in distances], c=colors[i], marker=("v" if iLabel==1 else "o"))
ax4.legend()
ax4.set_title('Mean Time for One Shot of a QEC Cycle')
ax4.set_xlabel('Distance')
ax4.set_ylabel('Time Elapsed')
ax4.grid(which='major')
In [ ]:
cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, len(distances)+1))
colors = [mplcolors.to_hex(c) for c in colors]

fig1, axs1 = plt.subplots(int(len(capacities)/2),2)
axs1=axs1.flatten()
fig1.set_size_inches( 2*fig4.get_size_inches()[0],fig4.get_size_inches()[1]*len(capacities)/1.4)

thresholds = []
thresholdIdxs = []
for j, (c, ax1) in enumerate(zip(capacities, axs1)):
    threshold = gate_improvements[0]
    thresholdIdx = 0
    for idx, improv in enumerate(gate_improvements[1:]):
        isThreshold = True
        for d1, d2 in zip(distances[:-1], distances[1:]):
            if data["LogicalErrorRates"]["Forwarding"][d1][j][idx+1]<data["LogicalErrorRates"]["Forwarding"][d2][j][idx+1]:
                isThreshold = False 
        if isThreshold:
            threshold = (gate_improvements[idx-1]+gate_improvements[idx])/2
            thresholdIdx = idx
            break
    thresholds.append(threshold)
    thresholdIdxs.append(thresholdIdx)

    for i,d in enumerate(distances):
        ax1.plot(gate_improvements, data["LogicalErrorRates"]["Forwarding"][d][j], linestyle='dashed', c=colors[i])
        ax1.scatter(gate_improvements, data["LogicalErrorRates"]["Forwarding"][d][j], c=colors[i], label=f"d={d}")
    
    ax1.vlines([threshold], ymin=[1e-7], ymax=[1e-0], colors=[colors[-1]], label="Threshold Point")
    ax1.legend()
    ax1.set_title(f'Threshold={threshold} for capacity={c}')
    ax1.set_xlabel('Gate Improvement')
    ax1.set_ylabel('Logical Error Rate')
    ax1.set_xlim(1e-2, 1000)
    ax1.set_ylim(1e-7, 1e-0)
    ax1.set_xticks(gate_improvements)
    ax1.loglog()

    ax1.grid(which='major')
In [ ]:
print(thresholds)
print(thresholdIdxs)
In [ ]:
thresholds[capacities.index(120)] = 7.5
thresholds[capacities.index(60)] = 7.5
thresholds[capacities.index(30)] = 7.5
thresholds[capacities.index(15)] = 7.5
thresholds[capacities.index(14)] = 7.5
thresholds[capacities.index(13)] = 7.5

for i, t in enumerate(thresholds):
    for j, g in enumerate(gate_improvements):
        if g>t:
            break 
    thresholdIdxs[i] = j
In [ ]:
for i in range(len(thresholds)):
    thresholds[i] = 7.5

for i, t in enumerate(thresholds):
    for j, g in enumerate(gate_improvements):
        if g>t:
            break 
    thresholdIdxs[i] = j
In [ ]:
cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, len(distances)+1))
colors = [mplcolors.to_hex(c) for c in colors]

fig1, axs1 = plt.subplots(int(len(capacities)/2),2)
axs1=axs1.flatten()
fig1.set_size_inches( 2*fig4.get_size_inches()[0],fig4.get_size_inches()[1]*len(capacities)/1.4)


for j, (c, ax1) in enumerate(zip(capacities, axs1)):
    threshold = thresholds[j]
    for i,d in enumerate(distances):
        ax1.plot(gate_improvements, data["LogicalErrorRates"]["Forwarding"][d][j], linestyle='dashed', c=colors[i])
        ax1.scatter(gate_improvements, data["LogicalErrorRates"]["Forwarding"][d][j], c=colors[i], label=f"d={d}")
    
    ax1.vlines([threshold], ymin=[1e-7], ymax=[1e-0], colors=[colors[-1]], label="Threshold Point")
    ax1.legend()
    ax1.set_title(f'Threshold={threshold} for capacity={c}')
    ax1.set_xlabel('Gate Improvement')
    ax1.set_ylabel('Logical Error Rate')
    ax1.set_xlim(1e-2, 1000)
    ax1.set_ylim(1e-7, 1e-0)
    ax1.set_xticks(gate_improvements)
    ax1.loglog()

    ax1.grid(which='major')
In [ ]:
def findMaxGi(_j, _d):
    return min(
        [
            u
            for u in range(len(gate_improvements))
            if data["LogicalErrorRates"]["Forwarding"][_d][_j][u] == 0
        ],
        default=len(gate_improvements),
    )


cs = np.concatenate(
    [
        np.concatenate(
            [
                [c] * len(gate_improvements[thresholdIdxs[j] : findMaxGi(j, d)])
                for d in distances
            ]
        )
        for j, c in enumerate(capacities)
    ]
)
ds = np.concatenate(
    [
        np.concatenate(
            [
                [d] * len(gate_improvements[thresholdIdxs[j] : findMaxGi(j, d)])
                for d in distances
            ]
        )
        for j in range(len(capacities))
    ]
)
xs = np.concatenate(
    [
        np.concatenate(
            [gate_improvements[thresholdIdxs[j] : findMaxGi(j, d)] for d in distances]
        )
        for j in range(len(capacities))
    ]
)
ys = np.concatenate(
    [
        np.concatenate(
            [
                [
                    data["LogicalErrorRates"]["Forwarding"][d][j][k + thresholdIdxs[j]]
                    for k in range(
                        len(gate_improvements[thresholdIdxs[j] : findMaxGi(j, d)])
                    )
                ]
                for d in distances
            ]
        )
        for j in range(len(capacities))
    ]
)

xs = np.log(xs)
ys = np.log(ys)


import sklearn.linear_model


_model2 = sklearn.linear_model.LinearRegression()


def _model2paramsfunc(_cs, _ds, _xs):
    _cs, _ds, _xs = np.array(_cs), np.array(_ds), np.array(_xs)
    return np.column_stack([ np.where(_cs==c, _xs, 0) for c in capacities[:-1]]+[np.where(_ds==d,_xs,0) for d in distances]+[np.where(_cs==c, 1, 0) for c in capacities[:-1]]+[np.where(_ds==d,1,0) for d in distances[:-1]])


def _model2rootfunc(_cs, _ds, _desired_y):
    _cs, _ds = np.array(_cs, dtype=int), np.array(_ds, dtype=int)
    _numer_inter = np.zeros_like(_cs, dtype=float)+p0_hat
    _denom_slope = np.zeros_like(_cs, dtype=float)
    _i = 0
    for _c in capacities[:-1]:
        _denom_slope+=ps_hat[_i]*np.where(_cs==_c, 1,0)
        _i+=1
    for _d in distances:
        _denom_slope+=ps_hat[_i]*np.where(_ds==_d,1,0)
        _i+=1
    for _c in capacities[:-1]:
        _numer_inter+=ps_hat[_i]*np.where(_cs==_c, 1,0)
        _i+=1
    for _d in distances[:-1]:
        _numer_inter+=ps_hat[_i]*np.where(_ds==_d,1,0)
        _i+=1

    return np.exp(
        np.divide(
            _desired_y - _numer_inter,
            _denom_slope,
        )
    )


_model2.fit(_model2paramsfunc(cs, ds, xs), ys)
p0_hat, ( ps_hat) = (
    _model2.intercept_,
    _model2.coef_,
)

print( p0_hat,  ps_hat)
ps_hat = ps_hat[2:]



_model2_errors_squared = (
    ys
    - _model2.predict(
        _model2paramsfunc(cs, ds, xs)
    )
) ** 2

print(np.sqrt(np.mean(_model2_errors_squared)))

x_bins = 100
y_bins = 100
_range = [[min(capacities), max(capacities)], [min(distances), max(distances)]]
counts, xedges, yedges = np.histogram2d(
    cs, ds, bins=[x_bins, y_bins], range=_range
)
values, _, _ = np.histogram2d(
    cs,
    ds,
    bins=[x_bins, y_bins],
    weights=_model2_errors_squared,
    range=_range,
)
counts = counts.T
values = values.T

with np.errstate(divide="ignore", invalid="ignore"):
    values = values / counts
    # scaled_to_error_values = (values/max(values.flatten()))*max(_model2_errors_squared)
    scaled_to_error_values = values
    scaled_to_error_values[counts == 0] = np.nan

x_centers = (xedges[:-1] + xedges[1:]) / 2
y_centers = (yedges[:-1] + yedges[1:]) / 2
X, Y = np.meshgrid(x_centers, y_centers)

points = np.column_stack([X.ravel(), Y.ravel()])
scaled_to_error_values_flat = scaled_to_error_values.ravel()

empty_bins = np.isnan(scaled_to_error_values_flat)


from scipy.interpolate import griddata
from scipy.spatial import cKDTree

if np.any(empty_bins):
    interpolated_values_flat = griddata(
        points[~empty_bins],
        scaled_to_error_values_flat[~empty_bins],
        points,
        method="cubic",
    )

empty_bins = np.isnan(interpolated_values_flat)
if np.any(np.isnan(interpolated_values_flat)):
    tree = cKDTree(points[~empty_bins])
    dist, idx = tree.query(points[np.isnan(interpolated_values_flat)])
    interpolated_values_flat[np.isnan(interpolated_values_flat)] = (
        interpolated_values_flat[~empty_bins][idx]
    )

interpolated_values = interpolated_values_flat.reshape(scaled_to_error_values.shape)
def _model2scalefunc(_cs, _ds):
    tree = cKDTree(points)
    _, idx = tree.query(list(zip(_cs, _ds)))
    return np.sqrt(interpolated_values_flat[idx])


print(np.mean(_model2scalefunc(cs, ds)))
In [ ]:
from scipy import stats
from scipy.optimize import minimize

cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, len(distances)+1))
colors = [mplcolors.to_hex(c) for c in colors]

fig1, axs1 = plt.subplots(int(len(capacities)/2),2)
axs1=axs1.flatten()
fig1.set_size_inches( 2*fig4.get_size_inches()[0],fig4.get_size_inches()[1]*len(capacities)/1.4)


fittings_dict = {c: {d: None for d in distances} for c in capacities}
for j, (c, ax1) in enumerate(zip(capacities, axs1)):
    for i,d in enumerate(distances):
        xs_for_cd = np.log(gate_improvements[thresholdIdxs[j]:findMaxGi(j,d)])
        if len(xs_for_cd)==0:
            continue
        ys = np.log(data["LogicalErrorRates"]["Forwarding"][d][j][thresholdIdxs[j]:findMaxGi(j,d)])
  
        res = stats.linregress(xs_for_cd, ys)
        alpha_hat, beta_hat, stderr_s_hat, stderr_i_hat = res.slope, res.intercept, res.stderr, res.intercept_stderr
        fittings_dict[c][d] = alpha_hat, beta_hat, (stderr_s_hat), (stderr_i_hat), np.sqrt(np.mean((ys-xs_for_cd*alpha_hat-beta_hat)**2))

        def readout_stat(_rys_drawn, _xs_proj):
            _res = stats.linregress(xs_for_cd, _rys_drawn)
            _alpha_hat, _beta_hat = _res.slope, _res.intercept
            return np.exp(_alpha_hat*(_xs_proj)+_beta_hat)
        
        def ry_star():
            return np.random.normal(loc=alpha_hat*xs_for_cd+beta_hat, scale=np.sqrt(stderr_s_hat*xs_for_cd+stderr_i_hat))

        xs_projection = np.log(np.array([1e-2, 1000]))
        readouts = np.array([readout_stat(ry_star(), xs_projection) for _ in range(10)]).T
        los, his = np.quantile(readouts, [0.025,0.975], axis=1)
  
        ax1.fill_between(np.exp(xs_projection), los, his, color=colors[i], interpolate=True,label=f"d={d}", linestyle="dashed", alpha=0.4)
        ax1.plot(np.exp(xs_projection), np.exp(alpha_hat*xs_projection+beta_hat), color=colors[i])
        ax1.scatter(gate_improvements, data["LogicalErrorRates"]["Forwarding"][d][j], c=colors[i], label=f"d={d}")
    
    ax1.vlines([thresholds[j]], ymin=[1e-7], ymax=[1e-0], colors=[colors[-1]], label="Threshold Point")
    ax1.legend()
    ax1.set_title(f'Threshold={threshold} for capacity={c}')
    ax1.set_xlabel('Gate Improvement')
    ax1.set_ylabel('Logical Error Rate')
    ax1.set_xlim(1e-2, 1000)
    ax1.set_ylim(1e-7, 1e-0)
    ax1.set_xticks(gate_improvements)
    ax1.loglog()

    ax1.grid(which='major')
In [ ]:
def rootfunc(_cs,_ds, _desired_y):
    _rs = []
    for c, d in zip(_cs, _ds):
        if not fittings_dict[c][d]:
            _rs.append(np.nan)
        else:
            _a_hat, _b_hat, _,_,_ = fittings_dict[c][d]
            _rs.append(np.exp((_desired_y-_b_hat)/_a_hat))
    return _rs

def _errorsSquared(_cs,_ds):
    _rs = []
    for c, d in zip(_cs, _ds):
        if not fittings_dict[c][d]:
            _rs.append(1.0)
        else:
            _, _, _, _, _mse = fittings_dict[c][d]
            _rs.append(_mse)
    return _rs

def _errorsFromFitting(_cs,_ds, _xs):
    _rs = []
    for c, d, x in zip(_cs, _ds,_xs):
        if not fittings_dict[c][d]:
            _rs.append(np.nan)
        else:
            _,_, _ss,_si,_ = fittings_dict[c][d]
            _rs.append(np.sqrt(_ss*x+_si))
    return _rs

def _predictFromFitting(_cs,_ds, _xs):
    _rs = []
    for c, d, x in zip(_cs, _ds,_xs):
        if not fittings_dict[c][d]:
            _rs.append(np.nan)
        else:
            _a_hat, _b_hat, _,_,_ = fittings_dict[c][d]
            _rs.append(_a_hat*x+_b_hat)
    return _rs
In [ ]:
cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, len(distances)+1))
colors = [mplcolors.to_hex(c) for c in colors]

desired_y = np.log(1e-4)

fig1, ax1 = plt.subplots()
for i, d in enumerate(distances[1:]):
    ax1.plot([k for k in range(len(capacities))], rootfunc(capacities, [d]*len(capacities), desired_y) , color=colors[i],label=f"d={d}")
    ax1.scatter([k for k in range(len(capacities))], rootfunc(capacities, [d]*len(capacities), desired_y), color=colors[i])
ax1.legend()
ax1.set_title(f'Improvement Needed to Survive a 10,000 Rounds')
ax1.set_xlabel('Capacities')
ax1.set_ylabel('Gate Improvement')

ax1.set_yticks(gate_improvements)
ax1.set_xticks([k for k in range(len(capacities))])
ax1.set_xticklabels(capacities)
ax1.set_yscale('log')

ax1.grid(which='major')
In [ ]:
cs = np.concatenate(
    [
        np.concatenate(
            [
                [c] * len(gate_improvements[thresholdIdxs[j] : findMaxGi(j, d)])
                for d in distances
            ]
        )
        for j, c in enumerate(capacities)
    ]
)
ds = np.concatenate(
    [
        np.concatenate(
            [
                [d] * len(gate_improvements[thresholdIdxs[j] : findMaxGi(j, d)])
                for d in distances
            ]
        )
        for j in range(len(capacities))
    ]
)
In [ ]:
capacities_proj = np.linspace(2,60,100)
distances_proj = np.linspace(3,9,100)
cs_proj = np.array([capacities_proj for _ in range(len(distances_proj))])
ds_proj = np.array([[d]*len(capacities_proj) for d in distances_proj])
desired_y = np.log(1e-4)
gate_improv_proj = np.log(rootfunc(cs, ds, desired_y))

_range = [[min(capacities), max(capacities)], [min(distances), max(distances)]]
_counts, cedges, dedges = np.histogram2d(
    cs, ds, bins=[100, 100], range=_range
)
_values, _, _ = np.histogram2d(
    cs,
    ds,
    bins=[100, 100],
    weights=gate_improv_proj,
    range=_range,
)
_counts = _counts.T
_values = _values.T

with np.errstate(divide="ignore", invalid="ignore"):
    scaled_to_improv_values = _values / _counts
    scaled_to_improv_values[_counts == 0] = np.nan

c_centers = (cedges[:-1] + cedges[1:]) / 2
d_centers = (dedges[:-1] + dedges[1:]) / 2
C, D = np.meshgrid(c_centers, d_centers)

_points = np.column_stack([C.ravel(), D.ravel()])
scaled_to_improv_values_flat = scaled_to_improv_values.ravel()

from scipy.interpolate import griddata
_empty_bins = np.isnan(scaled_to_improv_values_flat)
if np.any(_empty_bins):
    _interpolated_improv_values_flat = griddata(
        _points[~_empty_bins],
        scaled_to_improv_values_flat[~_empty_bins],
        _points,
        method="cubic",
    )

_interpolated_improv_values = _interpolated_improv_values_flat.reshape(scaled_to_improv_values.shape)
fig, axs = plt.subplots(1,2)
fig.set_size_inches(20,10)
ax=axs[0]
contour_filled = ax.contourf(C, D, np.exp(_interpolated_improv_values), cmap='viridis', norm=mplcolors.LogNorm())
ax.contour(C, D,np.exp( _interpolated_improv_values),  colors='black', norm=mplcolors.LogNorm())
cbar = fig.colorbar(contour_filled, ax=ax)
cbar.set_label('Gate Improvement')

ax.set_xlabel('Capacity')
ax.set_ylabel('Distance')

ax.set_yticks(distances)
ax.set_xlim(min(capacities),max(capacities[:-1]))
ax.set_ylim(min(distances),max(distances))

ax.set_title('Survive 10,000 Rounds')



x_bins = 100
y_bins = 100
_range = [[min(capacities), max(capacities)], [min(distances), max(distances)]]
counts, xedges, yedges = np.histogram2d(
    cs, ds, bins=[x_bins, y_bins], range=_range
)
values, _, _ = np.histogram2d(
    cs,
    ds,
    bins=[x_bins, y_bins],
    weights=np.exp(_errorsSquared(cs, ds)),
    range=_range,
)
counts = counts.T
values = values.T

with np.errstate(divide="ignore", invalid="ignore"):
    values = values / counts
    scaled_to_error_values = values
    scaled_to_error_values[counts == 0] = np.nan

x_centers = (xedges[:-1] + xedges[1:]) / 2
y_centers = (yedges[:-1] + yedges[1:]) / 2
X, Y = np.meshgrid(x_centers, y_centers)

points = np.column_stack([X.ravel(), Y.ravel()])
scaled_to_error_values_flat = scaled_to_error_values.ravel()

empty_bins = np.isnan(scaled_to_error_values_flat)


from scipy.interpolate import griddata
from scipy.spatial import cKDTree

if np.any(empty_bins):
    interpolated_values_flat = griddata(
        points[~empty_bins],
        scaled_to_error_values_flat[~empty_bins],
        points,
        method="cubic",
    )

empty_bins = np.isnan(interpolated_values_flat)
if np.any(np.isnan(interpolated_values_flat)):
    tree = cKDTree(points[~empty_bins])
    dist, idx = tree.query(points[np.isnan(interpolated_values_flat)])
    interpolated_values_flat[np.isnan(interpolated_values_flat)] = (
        interpolated_values_flat[~empty_bins][idx]
    )

interpolated_values = interpolated_values_flat.reshape(scaled_to_error_values.shape)




ax=axs[1]
contour_filled = ax.contourf( X, Y, np.exp(interpolated_values), cmap='viridis')
fig.colorbar(contour_filled, ax=ax, label='Squared Error')
ax.set_xlabel('Capacity')
ax.set_ylabel('Distance')

ax.set_yticks(distances)
ax.set_xlim(min(capacities),max(capacities[:-1]))
ax.set_ylim(min(distances),max(distances))
ax.set_title('Squared Error for Gate Improvement')
In [ ]:
_capacities = capacities[:-2]
cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, len(distances)))
colors = [mplcolors.to_hex(c) for c in colors]
desired_y = np.log(1e-4)
fig1, ax1 = plt.subplots()
for j, d in enumerate(distances[1:]):
    los, his = [], []
    for i,c in enumerate(_capacities):
        xs_for_cd = np.log(gate_improvements[thresholdIdxs[i]:findMaxGi(i,d)])
        if len(xs_for_cd)==0:
            los.append(0)
            his.append(0)
            continue
        def readout_stat(rys_drawn):
            _res = stats.linregress(xs_for_cd, rys_drawn)
            _alpha_hat, _beta_hat = _res.slope, _res.intercept
            return np.exp((desired_y-_beta_hat)/_alpha_hat)
        
        def ry_star():
            _cs, _ds = [c]*len(xs_for_cd), [d]*len(xs_for_cd)
            return np.random.normal(loc=_predictFromFitting(_cs, _ds, xs_for_cd), scale=_errorsFromFitting(_cs, _ds, xs_for_cd))

        readouts = np.array([readout_stat(ry_star()) for _ in range(10)])
        lo, hi = np.quantile(readouts, [0.025,0.975])
        los.append(lo)
        his.append(hi)   
    ax1.plot([k for k in range(len(_capacities))], rootfunc(_capacities, [d]*len(_capacities), desired_y) , color=colors[j],label=f"d={d}")
    ax1.scatter([k for k in range(len(_capacities))], rootfunc(_capacities, [d]*len(_capacities), desired_y), color=colors[j])
    ax1.hlines(los+his, xmin=[k-1/2 for k in range(len(_capacities))]*2, xmax=[k+1/2 for k in range(len(_capacities))]*2, colors=[colors[j]]*(len(los)+len(his)), alpha=0.8)
    ax1.vlines([k for k in range(len(_capacities))], ymin=los, ymax=his, colors=[colors[j]]*(len(los)+len(his)), alpha=0.8)
ax1.legend()
ax1.set_title(f'Improvement Needed to Survive 10,000 Rounds')
ax1.set_xlabel('Capacities')
ax1.set_ylabel('Gate Improvement')

ax1.set_xticks([k for k in range(len(_capacities))])
ax1.set_xticklabels(_capacities)
ax1.set_yscale('log')
ax1.grid(which='major')
In [ ]:
import scipy.stats

gi_offset = 0
_gate_improvements = gate_improvements[gi_offset:]
cmap = mpl.colormaps['plasma']
colors = cmap(np.linspace(0, 1, len(capacities)+1))
colors = [mplcolors.to_hex(c) for c in colors]

fig1, axs = plt.subplots(int(len(_gate_improvements)/2),2)
axs=axs.flatten()
fig1.set_size_inches( 2*fig4.get_size_inches()[0],fig4.get_size_inches()[1]*len(_gate_improvements)/1.4)
xs_projection = np.array([0, 50])

for j, (improv, ax) in enumerate(zip(_gate_improvements, axs)):
    for i, c in enumerate(capacities):
        if c==1 or c==30:
            logicalErrors = np.array([data["LogicalErrorRates"]["Forwarding"][d][i][j+gi_offset] for d in distances])
            ds = np.array(distances)
            xs, ys = ds[logicalErrors>0], logicalErrors[logicalErrors>0]
            ys = np.log(ys)

            def logPr(params, ys_given, eps = 1e-10):
                _alpha, _beta, _scale = params 
                _lik = stats.norm(loc=_alpha*(xs)+_beta, scale=_scale**2).pdf(ys_given)
                return np.log(_lik+eps)

            res = scipy.stats.linregress(xs, ys)
            alpha_hat, beta_hat, scale_s_hat, scale_i_hat = res.slope, res.intercept, res.stderr, res.intercept_stderr

            def readout_stat(rys_drawn, _xs):
                res = scipy.stats.linregress(xs, rys_drawn)
                _alpha_hat, _beta_hat = res.slope, res.intercept
                return np.exp(_alpha_hat*(_xs)+_beta_hat)
            
            def ry_star():
                return np.random.normal(loc=alpha_hat*xs+beta_hat, scale=np.sqrt(((ys-alpha_hat*xs-beta_hat)**2)))

            readouts = np.array([readout_stat(ry_star(), xs_projection) for _ in range(100)]).T
            los, his = np.quantile(readouts, [0.025,0.975], axis=1)
    
            ax.fill_between(xs_projection, los, his, color=colors[i], interpolate=True, label=f"c={c}", linestyle="dashed", alpha=0.4)


            ax.plot(xs_projection,
                np.exp(alpha_hat*xs_projection+beta_hat),
                linestyle='--',
                label=f"c={c}",
                c=colors[i])

for j, (improv, ax) in enumerate(zip(_gate_improvements, axs)):
    for i, c in enumerate(capacities):     
        if c==1 or c==30:
            logicalErrors = np.array([data["LogicalErrorRates"]["Forwarding"][d][i][j+gi_offset] for d in distances])  
            ds = np.array(distances)
            xs, ys = ds[logicalErrors>0], logicalErrors[logicalErrors>0] 
            ax.scatter(xs, ys, c=colors[i])
            
    ax.legend()
    ax.set_title(f'Projection for {improv}x Gate Improvement')
    ax.set_xlabel('Distance')
    ax.set_ylabel('Logical Error Rate')

    ax.set_ylim(1e-8, 1e2)
    ax.set_xlim(*xs_projection)
    ax.set_yscale('log')

    ax.grid(which='major')